The field of the invention is medical imaging and particularly, methods for reconstructing images from acquired image data.
In a computed tomography system, an x-ray source projects a fan-shaped beam which is collimated to lie within an x-y plane of a Cartesian coordinate system, termed the “image plane.” The x-ray beam passes through the object being imaged, such as a medical patient, and impinges upon an array of radiation detectors. The intensity of the transmitted radiation is dependent upon the attenuation of the x-ray beam by the object and each detector produces a separate electrical signal that is a measurement of the beam attenuation. The attenuation measurements from all the detectors are acquired separately to produce what is called the “transmission profile,” or “attenuation profile” or “projection.”
The source and detector array in a conventional CT system are rotated on a gantry within the imaging plane and around the object so that the angle at which the x-ray beam intersects the object constantly changes. The transmission profile from the detector array at a given angle is referred to as a “view” and a “scan” of the object comprises a set of views made at different angular orientations during one revolution of the x-ray source and detector. In a 2D scan, data is processed to construct an image that corresponds to a two dimensional slice taken through the object. The prevailing method for reconstructing an image from 2D data is referred to in the art as the filtered backprojection technique. This image reconstruction process converts the attenuation measurements acquired during a scan into integers called “CT numbers” or “Hounsfield units”, which are used to control the brightness of a corresponding pixel on a display.
The filtered backprojection image reconstruction method is the most common technique used to reconstruct CT images from acquired transmission profiles. As shown in FIG. 1 each acquired x-ray transmission profile 100 is backprojected onto the field of view (FOV) 102 by projecting each ray sum 104 in the profile 100 through the FOV 102 along the same ray path that produced the ray sum 104 as indicated by arrows 106. In projecting each ray sum 104 in the FOV 102 we have no a priori knowledge of the subject and the assumption is made that the x-ray attenuation in the FOV 102 is homogeneous and that the ray sum should be distributed equally in each pixel through which the ray path passes. For example, a ray path 108 is illustrated in FIG. 1 for a single ray sum 104 in one transmission profile 100 and it passes through N pixels in the FOV 102. The attenuation value, P, of this ray sum 104 is divided up equally between these N pixels:
      μ    n    =            (              P        ×        1            )        N  where μn is the attenuation value distributed to the nth pixel in a ray path having N pixels.
Clearly, the assumption that attenuation in the FOV 102 is homogeneous is not correct. However, as is well known in the art, if certain corrections are made to each transmission profile 100 and a sufficient number of profiles are acquired at a corresponding number of projection angles, the errors caused by this faulty assumption are minimized and image artifacts are suppressed. In a typical filtered backprojection method of image reconstruction, anywhere from 400 to 1000 views are typically required to adequately suppress image artifacts in a 2D CT image.
Another issue with x-ray CT is the x-ray dose to which the subject is exposed during the scan. To obtain a higher resolution and artifact free image it is necessary to obtain many views at a high enough x-ray beam intensity to reconstruct an image at the desired signal-to-noise ratio (SNR). The dose level may be reduced by decreasing the beam strength or reducing the number of acquired views, but either step also reduces the SNR of the reconstructed image.
Recent studies have underscored the importance of reducing the radiation dose to medical patients undergoing routine CT exams. For example, as reported by J. M. Peloquin, et al., in “Diagnostic Ionizing Radiation Exposure in a Population-Based Cohort of Patients with Inflammatory Bowel Disease,” American Journal of Gastroenterology 2008; 103:2015-2022, it was found that the radiation dose to patients with inflammatory bowel disease (IBD) was equivalent to, or greater than, the average radiation dose to the general population from naturally occurring sources, which is around 3.0 millisievert (mSv). Because as many as 1.4 million people in the United States and 2.2 million people in Europe are estimated to suffer from IBD, the development of imaging methods that reduce the radiation dose imparted to patients with IBD, while maintaining image quality, is important. In another recent study reported by J. E. Winslow, et al., in “Quantitative Assessment of Diagnostic Radiation Doses in Adult Blunt Trauma Patients,” Annals of Emergency Medicine 2008; 52:93-97, it was found that many trauma patients admitted to a Level I trauma center were exposed to clinically relevant levels of radiation exposure during their first 24 hours of care. This is because many trauma centers employ numerous x-ray imaging-based examinations when assessing trauma patients. Thus, it would be beneficial to develop a means for reducing the radiation dose imparted to trauma patients without affecting the quality of images obtained in the clinical setting.
Magnetic resonance imaging (MRI) uses the nuclear magnetic resonance (NMR) phenomenon to produce images. When a substance such as human tissue is subjected to a uniform magnetic field (polarizing field B0), the individual magnetic moments of the spins in the tissue attempt to align with this polarizing field, but precess about it in random order at their characteristic Larmor frequency. If the substance, or tissue, is subjected to a magnetic field (excitation field B1) which is in the x-y plane and which is near the Larmor frequency, the net aligned moment, MZ, may be rotated, or “tipped”, into the x-y plane to produce a net transverse magnetic moment Mxy. A signal is emitted by the excited spins, and after the excitation signal B1 is terminated, this signal may be received and processed to form an image.
When utilizing these signals to produce images, magnetic field gradients (Gx, Gy, and Gz) are employed. Typically, the region to be imaged is scanned by a sequence of measurement cycles in which these gradients vary according to the particular localization method being used. Each measurement is referred to in the art as a “view” and the number of views determines the quality of the image. The resulting set of received NMR signals, or views, or k-space samples, are digitized and processed to reconstruct the image using one of many well known reconstruction techniques. The total scan time is determined in part by the length of each measurement cycle, or “pulse sequence”, and in part by the number of measurement cycles, or views, that are acquired for an image. There are many clinical applications where total scan time for an image of prescribed resolution and SNR is a premium, and as a result, many improvements have been made with this objective in mind.
Projection reconstruction methods have been known since the inception of magnetic resonance imaging and this method is again being used as disclosed in U.S. Pat. No. 6,487,435. Rather than sampling k-space in a rectilinear, or Cartesian, scan pattern as is done in Fourier imaging and shown in FIG. 2A, projection reconstruction methods sample k-space with a series of views that sample radial lines extending outward from the center of k-space as shown in FIG. 2B. The number of views needed to sample k-space determines the length of the scan and if an insufficient number of views are acquired, streak artifacts are produced in the reconstructed image. The technique disclosed in U.S. Pat. No. 6,487,435 reduces such streaking by acquiring successive undersampled images with interleaved views and sharing peripheral k-space data between successive image frames.
Two example methods used to reconstruct images from an acquired set of projection views are described, for example, in U.S. Pat. No. 6,710,686. In MRI the most common method is to regrid the k-space samples from their locations on the radial sampling trajectories to a Cartesian grid. The image is then reconstructed by performing a 2D or 3D Fourier transformation of the regridded k-space samples. The second method for reconstructing an MR image is to transform the radial k-space projection views to Radon space by first Fourier transforming each projection view. An image is reconstructed from these signal projections by filtering and backprojecting them into the field of view (FOV). As is well known in the art, if the acquired signal projections are insufficient in number to satisfy the Nyquist sampling theorem, streak artifacts are produced in the reconstructed image.
Depending on the technique used, many MR scans currently used to produce medical images require many minutes to acquire the necessary data. The reduction of this scan time is an important consideration, since reduced scan time increases patient throughout, improves patient comfort, and improves image quality by reducing motion artifacts. Many different strategies have been developed to shorten the scan time.
One such strategy is referred to generally as “parallel imaging”. Parallel imaging techniques use spatial information from arrays of RF receiver coils to substitute for the encoding that would otherwise have to be obtained in a sequential fashion using RF pulses and field gradients (such as phase and frequency encoding). Each of the spatially independent receiver coils of the array carries certain spatial information and has a different sensitivity profile. This information is utilized in order to achieve a complete location encoding of the received MR signals by a combination of the simultaneously acquired data received from the separate coils. Specifically, parallel imaging techniques undersample k-space by reducing the number of acquired phase-encoded k-space sampling lines while keeping the maximal extent covered in k-space fixed. The combination of the separate MR signals produced by the separate receiver coils enables a reduction of the acquisition time required for an image (in comparison to conventional k-space data acquisition) by a factor that in the most favorable case equals the number of the receiver coils. Thus the use of multiple receiver coils acts to multiply imaging speed, without increasing gradient switching rates or RF power.
Two categories of such parallel imaging techniques that have been developed and applied to in vivo imaging are SENSE (SENSitivity Encoding) and SMASH (SiMultaneous Acquisition of Spatial Harmonics). With SENSE, the undersampled k-space data is first Fourier transformed to produce an aliased image from each coil, and then the aliased image signals are unfolded by a linear transformation of the superimposed pixel values. With SMASH, the omitted k-space lines are filled in or reconstructed prior to Fourier transformation, by constructing a weighted combination of neighboring lines acquired by the different receiver coils. SMASH requires that the spatial sensitivity of the coils be determined, and one way to do so is by “autocalibration” that entails the use of variable density k-space sampling.
The data acquisition methods are significantly different in the above exemplary imaging modalities. Namely, k-space is sampled to measure Fourier coefficients in MR data acquisitions while line integrals are measured in x-ray CT data acquisitions. Despite this, the challenge in image reconstruction for both modalities is common: reconstructing a high quality image from an undersampled data set.
According to standard image reconstruction theories, in order to reconstruct an image without aliasing artifacts, the sampling rate employed to acquire image data must satisfy the so-called Nyquist criterion, which is set forth in the Nyquist-Shannon sampling theorem. Moreover, in standard image reconstruction theories, no specific prior information about the image is needed. On the other hand, when some prior information about the desired or target image is available and appropriately incorporated into the image reconstruction procedure, an image can be accurately reconstructed even if the Nyquist criterion is violated. For example, if one knows a desired, target image is circularly symmetric and spatially uniform, only one view of parallel-beam projections (i.e., one projection view) is needed to accurately reconstruct the linear attenuation coefficient of the object. As another example, if one knows that a desired, target image consists of only a single point, then only two orthogonal projections that intersect at said point are needed to accurately reconstruct the image point. Thus, if prior information is known about the desired target image, such as if the desired target image is a set of sparsely distributed points, it can be reconstructed from a set of data that was acquired in a manner that does not satisfy the Nyquist criterion. Put more generally, knowledge about the sparsity of the desired target image can be employed to relax the Nyquist criterion; however, it is a highly nontrivial task to generalize these arguments to formulate a rigorous image reconstruction theory.
The Nyquist criterion serves as one of the paramount foundations of the field of information science. However, it also plays a pivotal role in modern medical imaging modalities such as magnetic resonance imaging (MRI) and x-ray computed tomography (CT). When the number of data samples acquired by an imaging system is less than the requirement imposed by the Nyquist criterion, artifacts appear in the reconstructed images. In general, such image artifacts include aliasing and streaking artifacts. In practice, the Nyquist criterion is often violated, whether intentionally or through unavoidable circumstances. For example, in order to shorten the data acquisition time in a time-resolved MR angiography study, undersampled projection reconstruction, or radial, acquisition methods are often intentionally introduced.
In contrast, undersampling is inevitable in four-dimensional cone beam CT (4D CBCT), such as when utilized in image-guided radiation therapy (IGRT). For example, in the case of IGRT, cone beam projection data are acquired over 10-15 respiratory cycles during a 60 second gantry rotation time. The acquired data is then retrospectively gated into 8-10 phases by synchronizing the respiratory signals with the data acquisition. After the respiratory gating, less than 100 cone beam projections are typically available to reconstruct images for each respiratory phase. Consequently, streaking artifacts are rampant in the reconstructed images for each respiratory phase. These undersampling artifacts pose a major challenge in 4D CBCT and limit the use of 4D CBCT in clinical practice.
Recently, a new image reconstruction method called highly constrained backprojection (HYPR) has been developed. As described in co-pending U.S. patent application Ser. No. 11/482,372, HYPR provides a method in which quality images can be reconstructed from far fewer projection signal profiles when a priori knowledge of the signal information is used in the reconstruction process. For example, signal information in an angiographic study may be known to include structures such as blood vessels. That being the case, when a backprojection path passes through these structures a more accurate distribution of a signal sample in each pixel can be achieved by weighting the distribution as a function of the known signal information at that pixel location. In HYPR, for a backprojection path having N pixels the highly constrained backprojection may be expressed as follows:
            S      n        =                  (                  P          ×                      C            n                          )                              ∑                      n            =            1                    N                ⁢                  C          n                      ,where Sn is the backprojected signal magnitude at a pixel n in an image frame being reconstructed, P is the signal sample value in the projection profile being backprojected, and Cn is the signal value of an a priori composite image at the nth pixel along the backprojection path. The composite image is reconstructed from data acquired during the scan, and may include that used to reconstruct the given image frame as well as other acquired image data that depicts the structures in the field of view. The numerator in the equation above, (P×Cn), weights each pixel using the corresponding signal value in the composite image and the denominator,
            ∑              n        =        1            N        ⁢          C      n        ,normalizes the value so that all backprojected signal samples reflect the projection sums for the image frame and are not multiplied by the sum of the composite image.
Also recently, a new mathematical framework for image reconstruction termed “compressed sensing” (CS) was formulated. In compressed sensing, only a small set of linear projections of a sparse image are required to reconstruct a quality image. The theory of CS is described in E. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information,” IEEE Transactions on Information Theory 2006; 52:489-509, and D. Donoho, “Compressed sensing,” IEEE Transactions on Information Theory 2006; 52:1289-1306, and is disclosed, for example, in U.S. patent application Ser. No. 11/199,675.
Although the mathematical framework of CS is elegant, the applicability of the reconstruction method to the field of medical imaging critically relies on the medical images being sparse. Unfortunately, a medical image is often not sparse in the standard pixel representation. Despite this, mathemtical transforms can be applied to a single image in order to sparsify the image. Such transforms are thus referred to as “sparsifying transforms”. More specifically, given a sparsifying transform, Ψ, CS image reconstruction is implemented by minimizing the following objective function:∥ΨI∥1,such that,AI=Y.
In the above objective function, I is a vector that represents the desired image, Y is a vector that represents the data acquired by the imaging system, A is a system matrix that describes the measurements, and
                            x                    1        =                  ∑                  i          =          1                N            ⁢                                x          i                              ,is the so-called L1-norm of an N-dimensional vector, x. Namely, the CS image reconstruction determines an image that minimizes the L1-norm of the sparsified image among all images that are consistent with the physical measurements, AX=Y.
The basic ideas in the compressed sensing image reconstruction theory can be summarized as follows. Instead of directly reconstructing a desired image in pixel representation, a sparsified version of the desired image is reconstructed. In the sparsified image, substantially fewer image pixels have significant image values; thus, it is possible to reconstruct the sparsified image from an undersampled data set. After the sparsified desired image is reconstructed, an “inverse sparsifying transform” is used to transform the sparsified image back to the desired image. In practice, there is no need to have an explicit form for the “inverse” sparsifying transform. In fact, only the sparsifying transform is needed in image reconstruction.